!>\brief
!!
!! Vlasov-Poisson equation.
!!
!! \n
!!
!! The electric field and potential are first diagnosed, and then the
!! RHS of the Vlasov equation is computed.
!!
!! We now consider the evaluation of the various terms appearing in
!! the Vlasov-Poisson equation. We assume that all the variables are
!! expressed in terms of the Lagrangian basis associated with the
!! Gauss-Lobatto quadrature points, and we assume that the coordinate
!! directions are orthogonal.
!!
!! The first term that we have to evaluate is
!! \f{displaymath}{
!!  \int_K f\underline{v}\cdot\nabla_x 
!!   \phi_{ \underline{i}_x\underline{i}_v }\,dx\,dv =
!!  \sum_{ \underline{l},\underline{m} }
!!    w_{\underline{l}} w_{\underline{m}} J_x J_v \,
!!    f(\underline{x}^G_{\underline{l}},
!!      \underline{v}^G_{\underline{m}})\,
!!    \underline{v}^G_{\underline{m}}\cdot\nabla_x
!! \phi_{ \underline{i}_x\underline{i}_v }
!!      (\underline{x}^G_{\underline{l}},
!!      \underline{v}^G_{\underline{m}}),
!! \f}
!! where \f$w\f$ denotes the quadrature weights,
!! \f$\underline{l},\underline{m}\f$ are multiindexes for space and
!! velocity, respectively, \f$J_x=|K_x|/|\hat{K}|\f$ (and the same for
!! \f$J_v\f$), and \f$\underline{x}^G,\underline{v}^G\f$ denote the
!! quadrature nodes in space and velocity, respectively. The
!! quadrature weights are readily obtained as
!! \f{displaymath}{
!!  w_{\underline{l}} = \prod_{i=1}^d w_{\underline{l}(i)}
!! \f}
!! and the coordinates of the quadrature nodes are obtained (up to a
!! suitable coordinate transformation) as
!! \f{displaymath}{
!!  \underline{x}^G_{\underline{l}} =
!!  \left[\begin{array}{c}
!!   x^G_{\underline{l}(1)} \\
!!   \vdots \\
!!   x^G_{\underline{l}(d)}
!!  \end{array}\right],
!! \f}
!! where it is understood that the one-dimensional weights and nodes
!! are considered when the subscript is a scalar. Less obvious is the
!! computation of the gradient of the test function. Let us first
!! notice that
!! \f{displaymath}{
!! \phi_{ \underline{i}_x\underline{i}_v }
!!      (\underline{x}, \underline{v}) =
!! \phi_{\underline{i}_x(1)}(x^1)\ldots
!! \phi_{\underline{i}_x(d)}(x^d)
!! \phi_{\underline{i}_v(1)}(v^1)\ldots
!! \phi_{\underline{i}_v(d)}(v^d), 
!! \qquad\underline{x}, \underline{v}\in\mathbb{R}^d,
!! \f}
!! so that
!! \f{displaymath}{
!! \frac{\partial}{\partial x^j}\phi_{ \underline{i}_x\underline{i}_v }
!!      (\underline{x}, \underline{v}) =
!! \phi_{\underline{i}_x(1)}(x^1)\ldots
!! \phi_{\underline{i}_x(j)}^\prime(x^j)\ldots
!! \phi_{\underline{i}_x(d)}(x^d)
!! \phi_{\underline{i}_v(1)}(v^1)\ldots
!! \phi_{\underline{i}_v(d)}(v^d), 
!! \qquad\underline{x}, \underline{v}\in\mathbb{R}^d.
!! \f}
!! Evaluating this relation at the quadrature points and including the
!! Jacobian \f$J_j=\frac{d x_j}{d\xi_j}\f$ yields
!! \f{displaymath}{
!! \frac{\partial}{\partial x^j}\phi_{ \underline{i}_x\underline{i}_v }
!!      (\underline{x}^G_{\underline{l}},
!!      \underline{v}^G_{\underline{m}}) =
!! \delta_{\underline{i}_x(1)\,\underline{l}(1)}\ldots
!! J_j^{-1}\phi_{\underline{i}_x(j)}^\prime(x^G_{\underline{l}(j)})\ldots
!! \delta_{\underline{i}_x(d)\,\underline{l}(d)}
!! \delta_{\underline{i}_v(1)\,\underline{m}(1)}\ldots
!! \delta_{\underline{i}_v(d)\,\underline{m}(d)}.
!! \f}
!! This implies that in the computation of the numerical integral only
!! \f$d\f$ terms must be included. To see this, let us introduce the
!! following notation for a generic multiindex \f$\underline{i}\f$:
!! \f{displaymath}{
!!  \underline{i}_{\backslash j} = \left[ \underline{i}(1)\,,
!!   \ldots, \underline{i}(j-1)\,,\,\underline{i}(j+1)\,,\ldots,
!!   \underline{i}(d) \right]
!! \f}
!! and
!! \f{displaymath}{
!!  \underline{i}_{\backslash j,k} = \left[ \underline{i}(1)\,,
!!   \ldots, \underline{i}(j-1)\,,k\,,\underline{i}(j+1)\,,\ldots,
!!   \underline{i}(d) \right].
!! \f}
!! Then, the only nonzero contributions are obtained when
!! \f{displaymath}{
!!  \underline{l}_{\backslash j} = 
!!  \underline{i}_{x_{\backslash j}},\qquad
!!  \underline{m} = \underline{i}_{v},
!! \f}
!! so that
!! \f{displaymath}{
!!  \int_K f\underline{v}\cdot\nabla_x 
!!   \phi_{ \underline{i}_x\underline{i}_v }\,dx\,dv =
!!    w_{\underline{i}_v} J_x J_v 
!!  \sum_{j=1}^d J_j^{-1}
!!    \underbrace{
!!    \left(\underline{v}^G_{\underline{i}_v}\right)_j}_{
!!     v^G_{\underline{i}_v(j)}}
!!  \sum_{k=1}^{p_k}
!!    w_{\underline{i}_{x_{\backslash j,k}}}
!!    f(\underline{x}^G_{\underline{i}_{x_{\backslash j,k}}},
!!      \underline{v}^G_{\underline{i}_v})\,
!! \phi_{\underline{i}_x(j)}^\prime(x^G_k),
!! \f}
!! Noting that
!! \f{displaymath}{
!!  w_{\underline{i}_{x_{\backslash j,k}}} = 
!!  w_{\underline{i}_{x_{\backslash j}}} w_k
!! \f}
!! and defining a differentiation matrix
!! \f{displaymath}{
!!  D_{ij}= w_j \phi_i^\prime(x^G_j),
!! \f}
!! we obtain
!! \f{displaymath}{
!!  \int_K f\underline{v}\cdot\nabla_x 
!!   \phi_{ \underline{i}_x\underline{i}_v }\,dx\,dv =
!!    w_{\underline{i}_v} J_x J_v 
!!  \sum_{j=1}^d J_j^{-1} v^G_{\underline{i}_v(j)}
!!  w_{\underline{i}_{x_{\backslash j}}}
!!  \sum_{k=1}^{p_k}
!!   D_{\underline{i}_x(j) k}\,
!!    f(\underline{x}^G_{\underline{i}_{x_{\backslash j,k}}},
!!      \underline{v}^G_{\underline{i}_v}),
!! \f}
!! where the second summation can be seen as the matrix multiplication
!! along the \f$j\f$-th spatial direction (contraction) of the
!! differentiation matrix (more precisely, of the
!! \f$\underline{i}_x(j)\f$-row of such matrix) and the local degrees
!! of freedom for \f$f\f$.
!!
!! Another term that must be computed is the boundary integral
!! \f{displaymath}{
!!  \int_{s_{v,x}}\underline{n}\cdot\widehat{\underline{v}f}
!!  \phi_{\underline{i}_x\underline{i}_v} \, d\sigma dv =
!!  \sum_{\underline{l}^s,\underline{m}}w_{\underline{l}^s}
!!  w_{\underline{m}} J^s_xJ_v\widehat{f}
!!  (\underline{x}^G_{\underline{l}^s},\underline{v}^G_{\underline{m}})
!!   \underline{v}^G_{\underline{m}}\cdot\underline{n}\,
!!   \phi_{\underline{i}_x\underline{i}_v}
!!   (\underline{x}^G_{\underline{l}^s},\underline{v}^G_{\underline{m}})
!! \f}
!! where \f$\underline{l}^s\f$ is a multiindex spanning the side
!! quadrature nodes. Notice that in fact there are two normals and two
!! test functions, one for each element connected to the side.
!! Clearly, this sum is zero unless \f$\underline{i}_x\in s_{v,x}\f$,
!! and in such case only one term must be considered. Moreover,
!! indicating by \f$n^s\f$ the normal direction of \f$s_{v,x}\f$, we
!! have
!! \f{displaymath}{
!!  w_{\underline{l}^s} = w_{\underline{l}_{\backslash n^s}},
!!  \qquad J^s = |s_{v,x}|, \qquad
!!  \underline{v}^G_{\underline{m}}\cdot\underline{n}= \pm
!!  v^G_{\underline{m}(n^s)}.
!! \f}
!! The sign of \f$\underline{v}^G_{\underline{m}}\cdot\underline{n}\f$
!! also defines the upwind flux \f$\widehat{f}\f$.
!!
!! \note For the use of \c p0 as initial guess for the solution of the
!! linear system, see comments in \c mod_poisson_pb.
!<----------------------------------------------------------------------
module mod_vp_ode

!-----------------------------------------------------------------------

 use mod_messages, only: &
   mod_messages_initialized, &
   error,   &
   warning, &
   info

 use mod_kinds, only: &
   mod_kinds_initialized, &
   wp

 !$ use omp_lib

 !$ use mod_omp_utils, only: &
 !$   mod_omp_utils_initialized, &
 !$   detailed_timing_omp, &
 !$   omput_push_key,    &
 !$   omput_pop_key,     &
 !$   omput_start_timer, &
 !$   omput_close_timer, &
 !$   omput_write_time

 use mod_state_vars, only: &
   mod_state_vars_initialized, &
   c_stv

 use mod_tps_phs_grid, only: &
   mod_tps_phs_grid_initialized, &
   t_tps_phs_grid, t_tps_grid

 use mod_tps_base, only: &
   mod_tps_base_initialized, &
   t_tps_base

 use mod_time_integrators, only: &
   mod_time_integrators_initialized, &
   c_ode, c_ods

 use mod_sparse, only: &
   mod_sparse_initialized, &
   ! sparse types
   t_col,       &
   t_tri,       &
   ! construction of new objects
   new_col,     &
   new_tri,     &
   ! convertions
   col2tri,     &
   tri2col,     &
   tri2col_skeleton, &
   tri2col_skeleton_part, &
   ! overloaded operators
   operator(+), &
   operator(*), &
   sum,         &
   transpose,   &
   matmul,      &
   ! error codes
   wrong_n,     &
   wrong_m,     &
   wrong_nz,    &
   wrong_dim,   &
   ! other functions
   nnz_col,     &
   nz_col,      &
   nz_col_i,    &
   get,         &
   set,         &
   diag,        &
   spdiag,      &
   ! deallocate
   clear

 use mod_octave_io_sparse, only: &
   mod_octave_io_sparse_initialized, &
   write_octave

 use mod_linsolver, only: &
   mod_linsolver_initialized, &
   c_linpb, c_itpb, c_mumpspb, c_umfpackpb, &
   gmres

 use mod_f_state, only: & 
   mod_f_state_initialized, &
   t_f_state, t_e_field

 use mod_poisson_pb, only: &
   mod_poisson_pb_constructor, &
   mod_poisson_pb_destructor,  &
   mod_poisson_pb_initialized, &
   compute_electric_field, compute_c11, p0

!-----------------------------------------------------------------------
 
 implicit none

!-----------------------------------------------------------------------

! Module interface

 public :: &
   mod_vp_ode_constructor, &
   mod_vp_ode_destructor,  &
   mod_vp_ode_initialized, &
   t_vp_ode, compute_diags, compute_electric_field, p0, &
   new_vp_ode, clear

 private

!-----------------------------------------------------------------------

! Module types and parameters

 !> Vlasov-Poisson ODE
 type, extends(c_ode) :: t_vp_ode
  type(t_tps_phs_grid), pointer :: grid    => null()
  type(t_tps_base),     pointer :: base(:) => null()
 contains
  procedure, pass(ode) :: rhs  => vp_rhs
 end type t_vp_ode
 
 interface clear
   module procedure clear_vp_ode
 end interface

! Module variables

 real(wp), allocatable :: m_mass(:,:,:), im_mass(:,:,:) !< mass matrix

 logical, protected ::               &
   mod_vp_ode_initialized = .false.
 character(len=*), parameter :: &
   this_mod_name = 'mod_vp_ode'

!-----------------------------------------------------------------------

contains

!-----------------------------------------------------------------------

 subroutine mod_vp_ode_constructor(grid,base,fu,alpha_rec)
  type(t_tps_phs_grid), intent(in), target :: grid
  type(t_tps_base), intent(in) :: base(2)
  integer, intent(in) :: fu
  logical, intent(in) :: alpha_rec
  character(len=*), parameter :: &
    this_sub_name = 'constructor'

  integer :: ie, iv, ix
  real(wp) :: xjw, vjw

   !Consistency checks ---------------------------
   if( (mod_kinds_initialized.eqv..false.) .or. &
    (mod_messages_initialized.eqv..false.) .or. &
!$      ( (detailed_timing_omp.eqv..true.).and. &
!$ (mod_omp_utils_initialized.eqv..false.) ) .or. &
  (mod_state_vars_initialized.eqv..false.) .or. &
(mod_tps_phs_grid_initialized.eqv..false.) .or. &
    (mod_tps_base_initialized.eqv..false.) .or. &
(mod_time_integrators_initialized.eqv..false.) .or. &
      (mod_sparse_initialized.eqv..false.) .or. &
(mod_octave_io_sparse_initialized.eqv..false.) .or. &
   (mod_linsolver_initialized.eqv..false.) .or. &
     (mod_f_state_initialized.eqv..false.) ) then  
     call error(this_sub_name,this_mod_name, &
                'Not all the required modules are initialized.')
   endif
   if(mod_vp_ode_initialized.eqv..true.) then
     call warning(this_sub_name,this_mod_name, &
                  'Module is already initialized.')
   endif
   !----------------------------------------------

   allocate( m_mass(base(1)%pkd,base(2)%pkd,grid%ne))
   allocate(im_mass(base(1)%pkd,base(2)%pkd,grid%ne))
   mass_do: do ie=1,grid%ne
     do iv=1,size(m_mass,2)
       vjw = grid%e(ie)%e1(2)%p%vol*base(2)%wgld(iv)
       do ix=1,size(m_mass,1)
         xjw = grid%e(ie)%e1(1)%p%vol*base(1)%wgld(ix)

          m_mass(ix,iv,ie) =     xjw*vjw
         im_mass(ix,iv,ie) = 1.0_wp/(xjw*vjw)
       enddo
     enddo
   enddo mass_do

   call mod_poisson_pb_constructor(grid,base(1),fu,alpha_rec)

   mod_vp_ode_initialized = .true.
 end subroutine mod_vp_ode_constructor

!-----------------------------------------------------------------------
 
 subroutine mod_vp_ode_destructor()
  character(len=*), parameter :: &
    this_sub_name = 'destructor'
   
   !Consistency checks ---------------------------
   if(mod_vp_ode_initialized.eqv..false.) then
     call error(this_sub_name,this_mod_name, &
                'This module is not initialized.')
   endif
   !----------------------------------------------

   call mod_poisson_pb_destructor()

   deallocate( m_mass)
   deallocate(im_mass)

   mod_vp_ode_initialized = .false.
 end subroutine mod_vp_ode_destructor

!-----------------------------------------------------------------------

 subroutine new_vp_ode(obj,grid,base)
  type(t_tps_phs_grid), intent(in), target :: grid
  type(t_tps_base),     intent(in), target :: base(2)
  type(t_vp_ode),       intent(out) :: obj

   obj%grid => grid
   obj%base => base

 end subroutine new_vp_ode
 
!-----------------------------------------------------------------------
 
 pure subroutine clear_vp_ode(obj)
  type(t_vp_ode), intent(out) :: obj

   nullify( obj%grid )
   nullify( obj%base )
 
 end subroutine clear_vp_ode
 
!-----------------------------------------------------------------------

 !> Compute the right-hand side of the Vlasov equation
 !!
 !! \note We assume that all the quadrature weights are defined on the
 !! interval \f$[0\,,1]\f$, so that the Jecobians are the element
 !! sizes.
 !!
 !! \note The mass matrix is not considered explicit: the differential
 !! operators already include the inverse of the diagonal mass matrix.
 !! This allows some speed-up of the code.
 subroutine vp_rhs(tnd,ode,t,uuu,ods,term)
  class(t_vp_ode), intent(in) :: ode
  real(wp),     intent(in)    :: t
  class(c_stv), intent(in)    :: uuu
  class(c_ods), intent(inout) :: ods
  class(c_stv), intent(inout) :: tnd
  integer,      intent(in), optional :: term(2)

  logical :: x_term, v_term
  integer :: ie, d, iv, ix, id, j, k1,k2,k3
  integer :: is, n, ixs, ivs, ix1, ix2, iv1, iv2, ie1, ie2, si1, si2
  integer, allocatable :: ixm(:), ivm(:), si(:)
  real(wp) :: &
    jsum, vjw, xjw, ivol1, ivol2, &
    egl,    & ! electric field at the quad. nodes (one component)
    vn, en, & ! normal velocity and electric field
    f1,f2,f_hat
  real(wp), allocatable :: vgl(:)
  !$ integer :: dotmap_work1(ode%base(1)%d), dotmap_work2(ode%base(2)%d)

  logical, parameter :: centered_xflux = .false.
  logical, parameter :: centered_vflux = .false.

   ! Check which terms are required
   x_term = .true.; v_term = .true. ! defaults
   if(present(term)) then
     if((term(1).gt.0).and.(term(2).gt.1)) then
       ! change the defaults
       if(term(1).le.term(2)/2) then
         v_term = .false.
       else
         x_term = .false.
       endif
     endif
   endif

   select type(uuu); type is(t_f_state)
     select type(tnd); type is(t_f_state)
       select type(ods); type is(t_e_field)

   ! Compute the electric field (store it in tnd)
   if(v_term) &
     call compute_electric_field(ods,uuu)

   tnd%f = 0.0_wp

   !$ if(detailed_timing_omp) then
   !$   call omput_push_key("vp_rhs_loops")
   !$   call omput_start_timer()
   !$ endif

   !$omp parallel &
   !$omp   private( d , ixm , ivm , vgl , si , ie , iv , id ,    &
   !$omp            ix , jsum , j , k1,k2,k3 , egl , is , n ,    &
   !$omp            ie1,ie2 , ivol1,ivol2 , vn , ixs , ix1,ix2 , &
   !$omp            xjw , f1,f2 , f_hat , si1,si2 , en , ivs ,   &
   !$omp            iv1,iv2 , vjw , dotmap_work1,dotmap_work2 )  &
   !$omp   shared( ode , x_term , v_term , uuu , ods , tnd )     &
   !$omp   default(none)

   d = ode%grid%d
   allocate(ixm(d),ivm(d),vgl(d),si(d))

   ! Element loop: there are essentially two contributions: advection
   ! in space and transport in velocity
   !$ if(detailed_timing_omp) then
   !$omp master
   !$   call omput_push_key("elem_do")
   !$   call omput_start_timer()
   !$omp end master
   !$ endif
   !$omp do schedule(static)
   elem_do: do ie=1,ode%grid%ne

     ie1 = ode%grid%e(ie)%e1(1)%p%o  ! space element
     si  = ode%grid%e(ie)%e1(2)%p%si ! sign indicator

     do iv=1,ode%base(2)%pkd ! Guass-Lobatto nodes
       ivm = ode%base(2)%mldx(iv,:) ! velocity multiindex
       if(x_term) then
        ! velocity at the quad. node
        do id=1,d
         vgl(id)=ode%grid%e(ie)%e1(2)%p%bw(id)*ode%base(2)%xgl(ivm(id))&
                  + ode%grid%e(ie)%e1(2)%p%box(1,id)
        enddo
        ! Include the bw prefactors
        vgl = vgl/ode%grid%e(ie)%e1(1)%p%bw
       endif
       do ix=1,ode%base(1)%pkd
         if(x_term) &
          ixm = ode%base(1)%mldx(ix,:) ! space multiindex

         jsum = 0.0_wp
         do j=1,d

           if(x_term) then
            ! f*v*grad_x(phi)
            call ode%base(1)%dotmap(k1,k2,k3 , ixm,j &
                                   !$ , dotmap_work1 &
                                   )
            jsum = jsum +                                              &
            ! The bw prefactor is moved outside the ix loop
            !vgl(j)/(ode%base%wgl(ixm(j))*ode%grid%e(ie)%e1(1)%p%bw(j))&
         vgl(j)/ode%base(1)%wgl(ixm(j))                                &
        * dot_product(ode%base(1)%wdphi(ixm(j),:),uuu%f(k1:k2:k3,iv,ie))
           endif

           if(v_term) then
            ! electric field at the quad. node
            egl = ods%et(j,si(j),ix,ie1)
            ! f*E*grad_v(phi)
            call ode%base(2)%dotmap(k1,k2,k3 , ivm,j &
                                   !$ , dotmap_work2 &
                                   )
            jsum = jsum -                                              &
            egl/(ode%base(2)%wgl(ivm(j))*ode%grid%e(ie)%e1(2)%p%bw(j)) &
        * dot_product(ode%base(2)%wdphi(ivm(j),:),uuu%f(ix,k1:k2:k3,ie))
           endif
         enddo
         ! quadrature weights are simplified with the mass matrix
         !jsum = vjw*xjw * jsum

         tnd%f(ix,iv,ie) = tnd%f(ix,iv,ie) + jsum
       enddo
     enddo

   enddo elem_do
   !$omp end do
   !$ if(detailed_timing_omp) then
   !$omp master
   !$   call omput_write_time()
   !$   call omput_close_timer()
   !$   call omput_pop_key()
   !$omp end master
   !$ endif

   ! x side loop: one contribution from advection in space
   if(x_term) then
    !$ if(detailed_timing_omp) then
    !$omp master
    !$   call omput_push_key("sidevx_do")
    !$   call omput_start_timer()
    !$omp end master
    !$ endif
    !$omp do schedule(static)
    sidevx_do: do is=1,ode%grid%nsvx

      n = ode%grid%sv_x(is)%s1%n ! normal space direction
      ie1 = ode%grid%sv_x(is)%ie(1)
      ie2 = ode%grid%sv_x(is)%ie(2)
      ivol1 = 1.0_wp/ode%grid%e(ie1)%e1(1)%p%vol
      ivol2 = 1.0_wp/ode%grid%e(ie2)%e1(1)%p%vol

      do iv=1,ode%base(2)%pkd
        ivm = ode%base(2)%mldx(iv,:) ! velocity multiindex
        vn  = ode%grid%sv_x(is)%e1%bw(n)*ode%base(2)%xgl(ivm(n)) &
             + ode%grid%sv_x(is)%e1%box(1,n)
        do ixs=1,ode%base(1)%spkd
          ! Geometry is obtained by left element (end element index)
          ix1 = ode%base(1)%s_dofs(2,n,ixs)
          ix2 = ode%base(1)%s_dofs(1,n,ixs)
          ixm = ode%base(1)%mldx(ix1,:) ! space multiindex
          xjw = ode%grid%sv_x(is)%s1%a/ode%base(1)%wgl(ixm(n))

          ! Numerical flux
          f1 = uuu%f( ix1,iv , ie1 )
          f2 = uuu%f( ix2,iv , ie2 )
          if(centered_xflux) then
            f_hat = 0.5_wp*(f1+f2)
          else
            if(vn.ge.0) then
              f_hat = f1
            else
              f_hat = f2
            endif
          endif

          ! Quadrature weights are simplified with the inverse mass
          ! matrix: this means that we only have to divide by the
          ! weight in the normal direction, stored in xjw, and consider
          ! the ratio surface/volume, using ivol.
          !tnd%f(ix1,iv,ie1) = tnd%f(ix1,iv,ie1) - vjw*xjw * vn*f_hat
          !tnd%f(ix2,iv,ie2) = tnd%f(ix2,iv,ie2) + vjw*xjw * vn*f_hat
          !$omp atomic
          tnd%f(ix1,iv,ie1) = tnd%f(ix1,iv,ie1) - xjw*ivol1 * vn*f_hat
          !$omp atomic
          tnd%f(ix2,iv,ie2) = tnd%f(ix2,iv,ie2) + xjw*ivol2 * vn*f_hat
        enddo
      enddo

    enddo sidevx_do
    !$omp end do
    !$ if(detailed_timing_omp) then
    !$omp master
    !$   call omput_write_time()
    !$   call omput_close_timer()
    !$   call omput_pop_key()
    !$omp end master
    !$ endif
   endif

   ! v side loop: one contribution from transport in velocty
   if(v_term) then
    !$ if(detailed_timing_omp) then
    !$omp master
    !$   call omput_push_key("sidexv_do")
    !$   call omput_start_timer()
    !$omp end master
    !$ endif
    !$omp do schedule(static)
    sidexv_do: do is=1,ode%grid%nsxv

      n = ode%grid%sx_v(is)%s1%n ! normal velocity direction
      ie1 = ode%grid%sx_v(is)%ie(1)
      ie2 = ode%grid%sx_v(is)%ie(2)
      ivol1 = 1.0_wp/ode%grid%e(ie1)%e1(2)%p%vol
      ivol2 = 1.0_wp/ode%grid%e(ie2)%e1(2)%p%vol
      si1 = ode%grid%e(ie1)%e1(2)%p%si(n)
      si2 = ode%grid%e(ie2)%e1(2)%p%si(n)

      do ix=1,ode%base(1)%pkd
        ! For the electric field we average the two values, which is
        ! significant if the sign indicators of the two elements are
        ! different.
        en = 0.5_wp * ( ods%et(n,si1,ix,ode%grid%sx_v(is)%e1%o) &
                      + ods%et(n,si2,ix,ode%grid%sx_v(is)%e1%o) )
        do ivs=1,ode%base(2)%spkd
          ! Geometry is obtained by left element (end element index)
          iv1 = ode%base(2)%s_dofs(2,n,ivs)
          iv2 = ode%base(2)%s_dofs(1,n,ivs)
          ivm = ode%base(2)%mldx(iv1,:) ! velocity multiindex
          vjw = ode%grid%sx_v(is)%s1%a/ode%base(2)%wgl(ivm(n))

          ! Numerical flux
          f1 = uuu%f( ix,iv1 , ie1 )
          f2 = uuu%f( ix,iv2 , ie2 )
          if(centered_vflux) then
            f_hat = 0.5_wp*(f1+f2)
          else
            if(en.le.0) then
              f_hat = f1
            else
              f_hat = f2
            endif
          endif

          !tnd%f(ix,iv1,ie1) = tnd%f(ix,iv1,ie1) + vjw*xjw * en*f_hat
          !tnd%f(ix,iv2,ie2) = tnd%f(ix,iv2,ie2) - vjw*xjw * en*f_hat
          !$omp atomic
          tnd%f(ix,iv1,ie1) = tnd%f(ix,iv1,ie1) + vjw*ivol1 * en*f_hat
          !$omp atomic
          tnd%f(ix,iv2,ie2) = tnd%f(ix,iv2,ie2) - vjw*ivol2 * en*f_hat
        enddo
      enddo

    enddo sidexv_do
    !$omp end do
    !$ if(detailed_timing_omp) then
    !$omp master
    !$   call omput_write_time()
    !$   call omput_close_timer()
    !$   call omput_pop_key()
    !$omp end master
    !$ endif
   endif

   deallocate(ixm,ivm,vgl,si)

   !$omp end parallel
   !$ if(detailed_timing_omp) then
   !$   call omput_write_time()
   !$   call omput_close_timer()
   !$   call omput_pop_key()
   !$ endif

       end select
     end select
   end select

 end subroutine vp_rhs

!-----------------------------------------------------------------------
 
 !> Compute the main diagnostics
 !!
 !! Concerning the electrostatic energy: \c e_l2norm is
 !! \f$\frac{1}{2}\int_{\Omega_x} \left( E_{x_1}^2 + \ldots +
 !! E_{x_d}^2 \right)dx\f$. More precisely, we use an array to
 !! separate the contributions for each direction: \f${\tt
 !! e\_l2norm(i)} = \frac{1}{2}\int_{\Omega_x} E_{x_i}^2 dx\f$. To
 !! obtain the total electrostatic energy, one simply takes
 !! <tt>sum(e_l2norm)</tt>. In principle, the same distinction could
 !! be made for the jumps of the electric potential, however they are
 !! typically very small.
 subroutine compute_diags(uuu , e_field , grid,base ,    &
   momentum,ekin,f_integr,f_l1norm,f_l2norm,f_min,f_max, &
   e_l2norm,p_jmnorm , ek)
  type(t_f_state), intent(in) :: uuu
  type(t_e_field), intent(inout), target :: e_field
  type(t_tps_phs_grid), intent(in) :: grid
  type(t_tps_base),     intent(in) :: base(2)
  real(wp), intent(out) :: momentum(:), ekin, f_integr, f_l1norm, &
    f_l2norm, f_min, f_max, e_l2norm(:), p_jmnorm
  real(wp), allocatable, intent(out) :: ek(:) ! spectral components
 
  integer :: ie, iv, is, id, ivm(grid%d), n, ie1, ie2, ix1, ix2, ixs
  real(wp) :: xvol_w(base(1)%pkd), xvol_wx(base(1)%pkdx)
  real(wp), dimension(base(2)%pkd) :: vvol_w, int_fx, int_afx, &
    int_f2x
  real(wp) :: vgl(grid%d), c11, sxa_w(base(1)%spkd)
  ! uuu_p is used to reshape the electric potential locally; see also
  ! section "12.6.6 Target Arguments" of the fortran handbook for the
  ! use of target dummy arguments associated with nontarget actual
  ! arguments.
  real(wp), pointer :: uuu_p(:,:)
  logical, parameter :: spectral_analysis = .true.
  character(len=*), parameter :: &
    this_sub_name = 'compute_diags'

   call compute_electric_field(e_field,uuu)

   ! Phase space integrals
   momentum = 0.0_wp
   ekin     = 0.0_wp
   f_integr = 0.0_wp
   f_l1norm = 0.0_wp
   f_l2norm = 0.0_wp
   f_min    =  huge(0.0_wp)
   f_max    = -huge(0.0_wp)
   elem_loop: do ie=1,grid%ne

     ! integration in x
     xvol_w = grid%e(ie)%e1(1)%p%vol * base(1)%wgld
     int_fx  = matmul( xvol_w ,     uuu%f(:,:,ie)    )
     int_afx = matmul( xvol_w , abs(uuu%f(:,:,ie))   )
     int_f2x = matmul( xvol_w ,     uuu%f(:,:,ie)**2 )

     ! integration in v
     vvol_w = grid%e(ie)%e1(2)%p%vol * base(2)%wgld
     do iv=1,base(2)%pkd ! Guass-Lobatto nodes
       ivm = base(2)%mldx(iv,:) ! velocity multiindex
       ! velocity at the quad. node
       do id=1,grid%d
         vgl(id) = grid%e(ie)%e1(2)%p%bw(id)*base(2)%xgl(ivm(id)) &
                  + grid%e(ie)%e1(2)%p%box(1,id)
         momentum(id) = momentum(id) + vvol_w(iv)*vgl(id)*int_fx(iv)
       enddo
       ekin = ekin + vvol_w(iv)*sum(vgl**2)*int_fx(iv) ! 1/2 added later
     enddo
     f_integr = f_integr + dot_product( vvol_w , int_fx  )
     f_l1norm = f_l1norm + dot_product( vvol_w , int_afx )
     f_l2norm = f_l2norm + dot_product( vvol_w , int_f2x ) ! sqrt later
     f_min    = min( f_min , minval(uuu%f(:,:,ie)) )
     f_max    = max( f_max , maxval(uuu%f(:,:,ie)) )

   enddo elem_loop
   ekin = 0.5_wp * ekin
   f_l2norm = sqrt(f_l2norm)

   ! space integrals: element integrals and side integrals
   e_l2norm = 0.0_wp
   x_elem_loop: do ie=1,grid%gx%ne
     do n=1,grid%d
       xvol_wx = grid%gx%e(ie)%vol * base(1)%wgldx(:,n)
       e_l2norm(n)=e_l2norm(n)+dot_product(xvol_wx,e_field%e(n,:,ie)**2)
     enddo
   enddo x_elem_loop
   ! Notice that more precisely we compute the electrostatic energy
   e_l2norm = 0.5_wp*e_l2norm

   uuu_p(1:base(1)%pkd,1:grid%gx%ne)         &
     => e_field%p%pl(1:base(1)%pkd*grid%gx%ne)
   p_jmnorm = 0.0_wp ! jump norm
   x_side_loop: do is=1,grid%gx%ns
     n = grid%gx%s(is)%n ! normal space direction
     ie1 = grid%gx%s(is)%ie(1)
     ie2 = grid%gx%s(is)%ie(2)
     c11 = compute_c11(is,grid%gx,base(1)%k)
     sxa_w = grid%gx%s(is)%a * base(1)%swgld(:,n)

     do ixs=1,base(1)%spkd
       ix1 = base(1)%s_dofs(2,n,ixs)
       ix2 = base(1)%s_dofs(1,n,ixs)
       p_jmnorm = p_jmnorm                                   &
         + sxa_w(ixs) * c11*(uuu_p(ix2,ie2)-uuu_p(ix1,ie1))**2
     enddo
   enddo x_side_loop
   ! Notice that, for consistency with e_l2norm, we introduce 1/2
   p_jmnorm = 0.5_wp*p_jmnorm

   ! Now the spectral components of the electric field: this is a
   ! preliminary implementation, and works only for 1D problems; also,
   ! one has to set manually the fundamental wave number.
   if(spectral_analysis) &
     call compute_spectral_coefficients(grid%gx,base(1),e_field%e,ek)

 contains

  pure subroutine compute_spectral_coefficients(gx,bx,e,ek)
   type(t_tps_grid), intent(in) :: gx
   type(t_tps_base), intent(in) :: bx
   real(wp), intent(in) :: e(:,:,:) ! electric field, nodal values
   real(wp), allocatable, intent(out) :: ek(:)
   
   integer, parameter :: nmodes = 4
   integer, parameter :: id = 1 ! only works in 1D
   !real(wp), parameter :: k0 = 0.5_wp ! Landau damping
   real(wp), parameter :: k0 = 0.3_wp ! bump-on-tail
   !real(wp), parameter :: k0 = 1.0_wp ! two-stream inst.
   real(wp), parameter :: & ! normalization coefficient
     scale = sqrt(k0/3.1415926535897932384626433832795028_wp)
   integer :: ie, in
   real(wp) :: k, sin_int(nmodes), cos_int(nmodes), xgx(bx%pkx)

    sin_int = 0.0_wp
    cos_int = 0.0_wp
    x_elem_loop: do ie=1,gx%ne
      ! nodes in physical space
      xgx = gx%e(ie)%bw(id)*bx%xglx + gx%e(ie)%box(1,id)
      do in=1,nmodes
        k = real(in,wp)*k0
        sin_int(in) = sin_int(in) + gx%e(ie)%bw(id) * dot_product( &
                                   bx%wglx , sin(k*xgx)*e(id,:,ie) )
        cos_int(in) = cos_int(in) + gx%e(ie)%bw(id) * dot_product( &
                                   bx%wglx , cos(k*xgx)*e(id,:,ie) )
      enddo
    enddo x_elem_loop

    allocate(ek(nmodes))
    ek = 0.5_wp * scale**2 * ( sin_int**2 + cos_int**2 )

  end subroutine compute_spectral_coefficients
 
 end subroutine compute_diags
 
!-----------------------------------------------------------------------

end module mod_vp_ode

